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Abstract. We investigate semi-discrete numerical schemes based on the stan- 
dard Galcrkin and lumped mass Galcrkin finite element methods for an initial- 
boundary value problem for homogeneous fractional diffusion problems with 
non-smooth initial data. We assume that f2 C R d , d = 1,2,3 is a convex 
polygonal (polyhedral) domain. We theoretically justify optimal order error 
estimates in L2- and i? 1 -norms for initial data in H~ s (£l), < s < 1. We 
confirm our theoretical findings with a number of numerical tests that include 
initial data v being a Dirac (5-function supported on a (d — l)-dimcnsional 
manifold. 

1. Introduction 

We consider the initial-boundary value problem for the fractional order parabolic 
differential equation for u(x, t): 

d?u(x,t) + Cu(x,t) = f(x,t), xintt T > t > 0, 

(1.1) u{x,t) = 0, xindn T>t>0, 

u(x, 0) = v(x), x in ft, 

where f2 C R d (d = 1,2,3) is a bounded convex polygonal domain with a boundary 
dfl, and £ is a symmetric, uniformly elliptic second-order differential operator. 
Integrating the second order derivatives by parts (once) gives rise to a bilinear form 
a(-, •) satisfying 

a{v,w) = (Cv,w) for all v e H 2 {VL),w e ^(fi), 

where (■,■) denotes the inner product in L,2(Cl). The form a(-,-) extends con- 
tinuously to Hq(Q) x Hq(Q) where it is symmetric and coercive and we take 
||«||ifi = a(u,u) 1 / 2 , for all u G Hq(Q). Similarly, C extends continuously to an 
operator from Hq(H) to i/~ 1 (f2) (the set of bounded linear functional on Hq(Q)) 

by 

(1.2) (Cu,v) = a(u,v) for all u,v e H^(Q). 



1991 Mathematics Subject Classification. 65M60, 65N30, 65N15. 

Key words and phrases, finite element method, fractional diffusion equation, error estimates, 
scmidiscretc discretization. 

The research of R. Lazarov and Z. Zhou was supported in parts by US NSF Grant DMS- 
1016525 and J. Pasciak has been supported by NSF Grant DMS-1216551. The work of all authors 
has been supported also by Award No. KUS-C1-016-04, made by King Abdullah University of 
Science and Technology (KAUST). 



2 



BANGTI JIN, RAYTCHO LAZAROV, JOSEPH PASCIAK, AND ZHI ZHOU 



Here (•,•) denotes duality pairing between H and H^fl). We assume that 

the coefficients of £ are smooth enough so that solutions v £ Hq(Q) satisfying 

a(v,4>) = (/,<£) for all 0e fl*(fi) 

with / e L 2 (fl) are in H 2 (Q). 

Here dfu (0 < a < 1) denotes the left-sided Caputo fractional derivative of order 
a with respect to t and it is defined by (cf. p. 91] or [TTJ p. 78]) 

i r*. . „ d 



9 ^ = ftMi (*-)" a ^)^ 

where T(-) is the Gamma function. Note that as the fractional order a tends to 
unity, the fractional derivative dfu converges to the canonical first-order derivative 



^ [§], and thus ( 1.1 1 reproduces the standard parabolic equation. The model ( 1.1 ) 



captures well the dynamics of subdiffusion processes in which the mean square 
variance grows slower than that in a Gaussian process PQ and has found a number 
of practical applications. A comprehensive survey on fractional order differential 
equations arising in viscoelasticity, dynamical systems in control theory, electrical 
circuits with fractance, generalized voltage divider, fractional-order multipoles in 
electromagnet ism, electrochemistry, and model of neurons is provided in [5]; see 
also HI]. 

The goal of this study is to develop, justify, and test a numerical technique for 



solving (1.1 1 with non-smooth initial data v G H~ S (Q), < s < 1, a important 
case in various applications and typical in related inverse problems; see e.g., [I], 
[l"2l Problem (4.12)] and [HE]. This includes the case of v being a delta- function 
supported on a (d — 1) -dimensional manifold in M. d , is particularly interesting from 
both theoretical and practical points of view. 



The weak form for problem (1.1) reads: find u{t) £ Hq(£1) such that 



(1.3) (d?u,x)+a(u, X ) = (f,x), Vye^(n), T>t>0, u(0) = 



The folowing two results are known, cf. [IS]: (1) for v € L 2 {£t) the problem ( 1.1 ) has 
a unique solution in C([0,T]; L 2 {VL) H C((0, T]; H 2 {Vl) n H^)) [12, Theorem 2.1]; 
(2) for / G L^O, T; L 2 (f2)), problem (Li} has a unique solution in L 2 (0, T; H 2 {n)n 
ffj(fl)) [H Theorem 2.2]. 



To introduce the semidiscrete FEM for problem (1.1) we follow standard no- 
tation in [2]. Let {Th} 0<h<1 be a family of regular partitions of the domain Q 
into d-simplexes, called finite elements, with h denoting the maximum diameter. 
Throughout, we assume that the triangulation Th is quasi-uniform, i.e., the diam- 
eter of the inscribed disk in the finite element r e Th is bounded from below by h, 
uniformly on Th- The approximation Uh will be sought in the finite element space 
Xh = Xh(p.) of continuous piecewise linear functions over Th'- 

Xh = {x & Hq(Q) : x is a linear function over t, Vt € Th} ■ 



The semidiscrete Galerkin FEM for problem (1.1) is: find Uh(t) e Xh such that 



(1.4) (a; Ma ) + aKx) = (/,x), 4 T>t >o, ^(0) = ^, 

where € X/j is an approximation of v. The choice of Vh will depend on the 
smoothness of v. For smooth data, v G H 2 (fl) n -Hq , we can choose u/j to 
be either the finite element interpolant or the Ritz projection onto Xh- In 
the case of non-smooth data, t> € £2(^)1 following Thomee [T3], we shall take 
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Vh = PhV, where Ph is the ^-orthogonal projection operator Ph : £2^) Xh, 
defined by (Ph&iX) — (0jX)j X S X h . In the intermediate case, v G Hq(Q), we 
can choose either Vh — PhV or Vh = The goal of this paper is to study 

the convergence rates of the semidiscrete Galerkin method ( |1.4[ ) for initial data 
v e H- S {n), < s < 1 when / = 0. 

The rest of the paper is organized as follows. In Section [2] we briefly review 
the regularity theory for problem . In Section [3] we motivate our study by 
considering a 1-D example with initial data being a 5-function. Then in Theorem 



3.1 we prove the main result: for < s < 1, the following error bound holds 

|K*) - u h (t)|| + h\\V(u(t) - u h (t))\\ < Ch 2 - s t- a £ h \\v\\^ s , 4 = | ln/i|. 

Further, in SectionHlwe show a similar result for the lumped mass Galerkin method. 
Finally, in Section p] we present numerical results for test problems with smooth, 
intermediate, non-smooth initial data and initial data that is a (5-function, all 
confirming our theoretical findings. 



2. Preliminaries 

For the existence and regularity of the solution to ( |1.1[ ), we need some notation 
and preliminary results. For s > — 1, we denote by H s (il) C i? _1 (r2) the Hilbcrt 
space induced by the norm 

oo 

(2-1) m^]Ta>,^) 2 

with {Xj}°Z 1 and {^jl^Li being respectively the Dirichlet eigenvalues and the 
L2-orthonormal eigenfunctions of L. As usual, we identify functions / in £2(0) 

with the functional F in ff-^fi) defined by (F, <fi) = for all (/> e H&(0,). 

1 

The set {<pj}JLij respectively, {AJ ^j}^!, forms an orthonormal basis in L 2 (fl), 
respectively, Thus |d|o = ||w|| = is the norm in £2(0) and — 

||«||H-i(n) is the norm in -ff -1 (f2). It is easy to check that |w|i = a(v, v) 2 is also the 
norm in _ffg(J7). Note that {H s (il)}, s > — 1 form a Hilbert scale of interpolation 
spaces. Motivated by this, we denote || • \\h° to be the norm on the interpolation 
scale between H^fl) and L 2 (£l) when s is in [0, 1] and || • to be the norm on the 
interpolation scale between L2^l) and if _1 (fi) when s is in [—1,0]. Thus, || • \\h" 
and I • | s provide equivalent norms for s € [—1, 1]. 

We further assume that the coefficients of the elliptic operator C are sufficiently 
smooth and the polygonal domain f2 is convex, so that | v| 2 = ||£w|| is equivalent to 
the norm in H 2 {n) n H£(Q) (cf. the proof of Lemma 3.1 of [T^). 

Now we introduce the operator E{t) by 

00 00 

(2.2) E(t)v = J2E<xA-^t a ) {v,<Pi)<Pi, ^ere E a . (z) = T(k l +8) - 
j=l k=o [ + 1 ' 

Here E a ^p(z) is the Mittag-Leffler function defined for z € C [§]. The operator E(t) 



gives a representation of the solution u of (1.1 1 with a homogeneous right hand 
side, so that for f(x,t) = we have u(t) = E(t)v. This representation follows from 
cigenfunction expansion [T2]. Further, we introduce the operator E(t) defined for 
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X <E L 2 (f2) as 

oo 

(2.3) E(t)x = ^t a - 1 E a<a (-\ j t a ) ( X ,<Pj)<Pj- 

j=o 

The operators E(t) and E(t) together give the following representation of the so- 
lution of ( [Llj ): 

(2.4) u(t)=E(t)v+ [ E(t- s)f{s)ds. 



Motivated by [4j [12], we will study the convergence of semidiscrete Galerkin 
methods for problem ( |1.1[ ) with very weak initial data, i.e., v € H~ S (Q), < s < 1. 
Then the following question arises naturally: in what sense should we understand 
the solution for such weak data? Obviously, for any t > the function u(t) = E(t)v 
satisfies equation (1.1 1. Moreover, by dominated convergence we have 

oo l 

]ha + \E(t)v-v\. s = ( t lim £(J5 ail (-AjO " l^'K^) 2 )* = 

3=1 

provided that v G H~ S (Q). Here (v,Lpj) = (v,ipj)H-=,H 3 is we U defined since 
ifj € Hq(£1). Therefore, the function u(t) — E(t)v satisfies (1.1) and for £ — J- it 
converges to v in H~ s -noim. That is, it is a weak solution to (1.1); see also [4j 
Proposition 2.1]. 

For the solution of the homogeneous equation ( |1.1| , which is the object of our 
study, we have the following stability and smoothing estimates. 



Theorem 2.1. Let u{t) = E(t)v be the solution to problem (1.1 ) with / = 0. Then 
for t > we have the the following estimates: 

(a) for £ = 0, < q < p < 2 and for £ = I, < p < q < 2 and q < p + 2: 

(2.5) mfu{t)\ p <Ct'^ + ^\v\ q , 

(b) for < s < 1 and < p + s < 2 

(2.6) \d?u(t)\- s <Ci->|_ s , and \u(t)\ p < Ct-^ a \v\_ s . 

Proof. Part (a) can be found in J2J Theorem 2.1] and [B] Theorem 2.1]. Hence we 
only show part (b). Note that for t > 0, 

oo oo 

mil <Y,y\E a A~^n\ 2 \(v,<Pi)\ 2 < (1 + x. ta) 2 \^^)\ 2 

3=0 3=0 \ 3 ' 

( i + x 3 ta ) 

oo 

< Ct-^ a Y,^\{vA 3 )\ 2 = Ct~^ a \vt s , 

3=0 

which proves the second inequality of case (b). The first estimate follows similarly 
by noticing the identity d?E aA (-\t a ) = -XE aA (-Xt a ) [§]. □ 

We shall need some properties of the L 2 -projection onto X^. 
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Lemma 2.1. Assume that the mesh is quasi-uniform. Then for s G [0, 1], 

- P h )w\\ H s < Ch 2 - s \\w\\ H 2, for all w G H 2 (fl) n H^(Cl), 

and 

||(7- P h )w\\ H . < C7i 1_s |HlH l i for allw e Hl(VL). 
In addition, Ph is stable on H s (fl) for s £ [—1,0]. 

Proof. Since the mesh is quasi-uniform, the /^-projection operator Ph is stable in 
i?o(f2) l 2J. This immediately implies its stability in Thus, stability on 

H~ s (fl) follows from this, the trivial stability of Ph on £2(0) an d interpolation. 

Let J/j be the finite element interpolation operator and Ch be the Clement or 
Scott-Zhang interpolation operator. It follows from the stability of Ph in £2^) 
and H$(Q) that 

\\{I - P h )w\\ L2 < \\(I-I h )w\\ L2 <Ch 2 \\w\\ H 2, for all we H 2 (Q) n H^(Q), 

- P h )w\\ H i < C\\(I - I h )w\\ H i <Ch\\w\\ H 2, for all w G H 2 (n) n H^(Q), 

- Ph)w\\ L2 < - C h )w\\ L2 <Ch\\w\\ H i, for all ^Gi^), 
P h )w\\m < C\\w\\ H i, for all w G H^(Q). 

The inequalities of the lemma follow by interpolation. □ 



Remark 2.2. All the norms appearing in Lemma \2.1\ can be replaced by their 
corresponding equivalent dotted norms. 

3. GALERKIN FINITE ELEMENT METHOD 

To motivate our study we shall first consider the 1-D case, i.e., Cu = —u", and 
take initial data the Dirac <5-function at x = \, (S,v) — v(^). It is well known that 

Hq (0,1) embeds continuously into Co(0, 1), hence the (5-function is a bounded 

l_|_ e 

linear functional on the space °(fi), i.e., 5 G H-2-\tt). 

In Tables[l]and[2] we show the error and the convergence rates for the semidiscrete 
Galerkin FEM and semidiscrete lumped mass FEM (cf. Section 4) for initial data 
v being a Dirac (5-function at x = \. The results suggest an 0(/i5) and 0(h2~) 
convergence rate for the H 1 - and LVnorm of the error, respectively. Below wc 
prove that up to a factor | lnh\ for fixed t > 0, the convergence rate is of the order 
reported in Tables [l] and [2] In Table [3] we show the results for the case that the 
(5-function is supported at a grid point. In this case the standard Galerkin method 
converges at the expected rate in -H^-norm, while the convergence rate in the L2- 
norm is 0(h 2 ). This is attributed to the fact that in 1-D the solution with the 
(5-function as the initial data is smooth from both sides of the support point and 
the finite element spaces have good approximation property. 

The numerical results in Tables [T]-[3] motivate our study on the convergence rates 
of the semidiscrete Galerkin and lumped mass schemes for initial data v G H~ S (VL) : 
< s < 1. 



Theorem 3.1. Letu anduh be the solutions of (1.1) and the semidiscrete Galerkin 



finite element method (1.4) with Vh — PhV, respectively. Then there is a constant 



C > such that for < s <1 

(3.1) ||« fc (t) - u(t)|| + h\\V(u h (t) - u(t))\\ < Ch 2 - S 4 t~ a \v\. 
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Table 1. Standard FEM with initial data for h = l/(2 k + l), 
a = 0.5. 



time 


k 


3 


4 


5 


6 


7 


ratio 


rate 


t = 0.005 


L 2 " n orm 


3.95e-2 


1.59e-2 


6.00e-3 


2.19e-3 


7.89e-4 


ps 2.75 


O(0) 


if -norm 


1.21e0 


8.99e-l 


6.52e-l 


4.66e-l 


3.33e-l 


ps 1.40 


0(h*) 


t = 0.01 


Z/2-norm 


2.85e-2 


1.13e-2 


4.26e-3 


1.55e-3 


5.58e-4 


ps 2.77 


O(0) 


i/^-norm 


8.66e-l 


6.39e-l 


4.62e-l 


3.31c-l 


2.35e-l 


pa 1.40 


O(hi) 


t = 1 


L2-norm 


3.04c-3 


1.17e-3 


4.34e-4 


1.57e-4 


5.61e-5 


ps 2.79 


0(hl) 


H -norm 


8.91e-2 


6.49e-2 


4.66e-2 


3.32e-2 


2.36e-2 


ss 1.41 


O(hi) 


Table 2. Lumped mass FEM with initial data <J(|), h = l/2 k a = 0.5. 


time 


k 


3 


4 


5 


6 


7 


ratio 


rate 


t = 0.005 


Z/2 _n orm 


7.24e-2 


2.66e-2 


9.54e-3 


3.40e-3 


1.21e-3 


ps 2.79 


O(ftt) 


H -norm 


1.51c0 


1.07e0 


7.60e-l 


5.40e-l 


3.81e-l 


PS 1.41 


0(h*) 


t = 0.01 


L2-norm 


5.20e-2 


1.89e-2 


6.77e-3 


2.40e-3 


8.54e-4 


ps 2.79 


0(0) 


iJ 1 -norm 


1.07e0 


7.59e-l 


5.37e-l 


3.80e-l 


2.70e-l 


ps 1.41 


0(hi) 


t = 1 


Z/2-norm 


5.47e-3 


1.93e-3 


6.84e-4 


2.42e-4 


8.56e-5 


ps 2.79 


O(0) 


if 1 -norm 


1.07e-l 


7.58e-2 


5.37e-2 


3.80e-2 


2.70e-2 


ps 1.41 


0(h*) 



Table 3. Standard semidiscrete FEM with initial data h — 

l/2 k , a = 0.5. 



Time 




3 


4 


5 


6 


7 


ratio 


rate 


t = 0.005 


Z/2 _n orm 


5.13e-3 


1.28e-3 


3.21e-4 


8.03e-5 


2.01e-5 


ps 3.99 


0(h 2 ) 


H -norm 


4.29e-l 


3.09e-l 


2.21e-l 


1.56e-l 


l.lle-1 


ps 1.41 


0{hi) 


t = 0.01 


I/2 _n orm 


3.07e-3 


7.70e-4 


1.93e-4 


4.82e-5 


1.21e-5 


ps 3.98 


0(h 2 ) 


-norm 


3.04e-l 


2.19e-l 


1.56e-l 


l.lle-1 


7.87e-2 


ps 1.41 


0(hi) 


t = 1 


Z/2 _n orm 


1.44e-5 


2.64e-6 


6.66e-7 


1.69e-7 


4.30e-8 


ps 3.94 


0(h 2 ) 


H -norm 


3.15e-2 


2.23e-2 


1.58e-2 


l.lle-2 


7.81e-3 


ps 1.41 


0{hi) 



Remark 3.2. Note that for any fixed e there is a C e > such that |5|_i_ e < C e . 
Thus, modulo the factor £h = \ hih\, the theorem confirms the computational results 
of Tableu^ namely convergence in the L2~norm with a rate 0(hi) and in H 1 -norm 
with a rate OQi*). 

Proof. We shall need the following auxiliary problem: find u h (t) € ZZg(fi), s.t. 
(3.2) (d?u h (t),x) + a(u h (t), X ) = (f(t),x), V X e H^Q), t > 0, u h (0) = P h v. 
We note that the initial data u h (0) = P^v E Hq(CI) is smooth. 
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Now we consider the semidiscrete Galerkin method for problem (3.2), i.e., equa- 
th v h - 

u h (t) - if 



tion (1.4 1 with Vh = P^v. By Theorem 3.2 of [6] we have 
(3.3) \\u h (t)-u h (t)\\ + h\\V(u h (t)-u h (t))\\<Ch 2 £ h t- a \\P h v\\. 

for < s < 1, and the 



Now, using the inverse inequality ||P/ l 'u|] < Ch s ||P/ l u| 
stability of Ph in H~ s (n) (cf. Lemma 2.1 1, we get 

(3.4) \\u h (t)-u h (t)\\+h\\V(u h (t)-u h (t))\\<Ch 2 - s £ h t- a \\v\\. s . 

Now we estimate u(t) — u h (t) — E{t)(v — Phv). To this end, let {v n } c L 2 {£1) 
be a sequence converging to v in H~ s (tt). Noting that the operators Ph and E{t) 
are self-adjoint in (•, •) and using the smoothing property (2.5| of E{t) with t = 0, 
q = and p = 2, we obtain for any <j> 6 L^(fl) 

\{E{t){I - P h )v n ,4>)\ = \(v n ,(I-P h )E(t)<P)\ < \v n \- s \(I - P h )E(t)cj>\ s 

< Ch 2 - s \v n \^ s \E{t)^\ 2 < Ch % -H- a \v n \_ s \\<j>l 

Taking the limit as n tends to infinity gives 

II (A fc^ll \{E{t){I - P h )v ^)\ 

3.5 \\u{t)-u{t)\\= sup jr-jj 

4>£L 2 (n) \\n 



< ch 2 - s r a \v\ 



Then by the triangle inequality we arrive at the L 2 -estimate in (3.1). 

Next, for the gradient term || V(u(t)— u h (t))\\, we observe that for any G 7J 1 (r2), 
by the coercivity of a(-, •), we have 

C \\V(E{t)(I - P h )v n )\\ 2 < a(E(t)(I - P h )v n , E(t)(I - P h )v n ) 
(3.6) < ^ a(E{t){I ~ P h )v n ,4>) 2 

~ 4>eHl(Q.) a,{^,4>) 

Meanwhile we have 

\a(E{t)(I-P h )v n ,c(>)\ = \((I-P h )v n ,E(t)C<t>)\ = \(v n ,(I-P h )E(t)£4>)\ 

< C\vn\- a \(I - P h )E(t)£<P\ s < Ch x - s \v n \- s \E{t)C<l>\i 

< Ch^H-^v^C^x < Ch l - s t- a \v n \- a \(j>\i- 



Passing to the limit as n tends to infinity and combining with (3.6) gives 
(3.7) ||V(«(*) - u fc (t))|| < Ctf-H-^vls. 



Thus, ( |3.5| and (3.7) lead to the following estimate for < s < 1: 
(3.8) ||«(t) - u fc (t)|| + h\\V(u(t) - u h (t))\\ < Ch 2 - S r a \v\. 



Finally, (3.4), (3.8), and the triangle inequality give the desired estimate (3.1) and 
this completes the proof. □ 



4. Lumped mass method 

In this section, we consider the lumped mass FEM in planar domains (see, e.g. 
[21 Chapter 15, pp. 239-244]). An important feature of the lumped mass method is 
that when representing the solution tih in the nodal basis functions, the mass matrix 
is diagonal. This leads to a simplified computational procedure. For completeness 
we shall briefly describe this approximation. Let zj, j = 1, . . . , d + 1 be the vertices 
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of the (i-simplex r G Th ■ Consider the following quadrature formula and the induced 
inner product in X^: 

QrAf) = E ~ / f dx > ( w ^)h = E QrA^x) 

3=1 Jt rdTh 

Then lumped mass finite clement method is: find Uh{t) G X/, such that 
(4.1) (d?u hlX )h+a(u h , X ) = if,X) V X G X h , t > 0, u h (0) = P h v. 

To analyze this scheme we shall need the concept of symmetric meshes. Given 
a vertex z G Th, the patch II z consists of all finite elements having z as a vertex. 
A mesh Th is said to be symmetric at the vertex z, if x G H z implies 2z — x G II Z , 
and 7ft is symmetric if it is symmetric at every interior vertex. 

In [6j Theorem 4.2] it was shown that if the mesh is symmetric, then the lumped 



mass scheme (4.1) for / = has an almost optimal convergence rate in L 2 -norm 
for nonsmooth data v G L,2(p,)- 

Now we prove the main result concerning the lumped mass method: 



Theorem 4.1. Let u(t) and Uh(t) be the solutions of the problems (1.1) and (4.1), 

respectively. Then for t > the following error estimate is valid: 

(4.2) ||fi h (t) - u(t)|| + ||V(« h (t) - u(t))|| < Ch 1 -'t h t- a \\v\\- a , < s < 1. 
Moreover, if the mesh is symmetric then 

(4.3) ||u fc (t) - u(t)|| < Ch 2 - s e h t- a \\v\\- s , < s < 1. 

Proof. We split the error into Uh(t) — u(t) = Uh(t) — u h (t) + u h (t) — u(t), where 
u h (t) — u(t) was estimated in ( |3.8| ). The term Uh(t)~u h (t) is the error of the lumped 
mass method for the auxiliary problem ( |3.2[ ). Since the initial data P^v G -L 2 (f2), 
we can apply known estimates on Uh(t) — vr(t) [SJ Theorem 4.2]. Namely, 

(a) If the mesh is globally quasiuniform, then 

||ti fc (f) - u h (t)\\ + h\\W(u h (t) - u h (t))\\ < Chr a £ h \\P h v\\; 

(b) If the mesh is symmetric, then 

\\u h (t)-u h (t)\\ < ChH- a £ h \\P h v\\. 
These two estimates, the inequality ||P/iu|| < C/i _s ||w||_ s , < s < 1, and estimate 



(3.4) give the desired result. This completes the proof of the theorem. □ 



Remark 4.2. The H 1 -estimate is almost optimal for any quasi-uniform meshes, 
while the L^- estimate is almost optimal for symmetric meshes. For the standard 
parabolic equation with initial data v G L 2 (Q), it was shown in [3] that the lumped 
mass scheme can achieve at most an 0(h*) convergence order in L 2 -norm for some 
nonsymmetric meshes. This rate is expected to hold for fractional order differential 
equations as well. 

5. Numerical results 
Here we present numerical results in 2-D to verify the error estimates derived 



herein and [B]. The 2-D problem (1.1) is on the unit square £1 = (0, l) 2 with 



C = — A. We perform numerical tests on four different examples: 
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(a) Smooth initial data: v(x, y) — x(l — x)y(l — y)\ in this case the initial data 
v is in H 2 (Vl) n H^fl), and the exact solution u(x,t) can be represented 
by a rapidly converging Fourier series: 

u(x,t) = X] X! 3 U 3 m e E a,i(-^n, m t a ) sin(mra) sm(rmry), 

n—l m—1 

where A„. m = (n 2 + to 2 )7t 2 , and q = 4sin 2 (^7r/2) — i7rsin(Z7r), I = m,n. 

(b) Initial data in Hq(CI) (case of intermediate smoothness): 

v(x) = (x- \){x - l)(y - \){y - l)X[i,i] x [i,i], 

where X[i i]x[A 1] ^ s the characteristic function of [g, 1] x [|, 1]. 

(c) Nonsmooth initial data: v(x) = X[i §]x[A §]■ 

(d) Very weak data: v = Sr with T being the boundary of the square [4, |] x 
[|, |] with (Sr, 4>) = J r <fr(s)ds. One may view (v, x) for x^ftC i?3+ e (f2) 
as duality pairing between the spaces H~^~ € (£l) and H^ +e (Vt) for any e > 
so that <5r G -ff~ 5_£ (0). Indeed, it follows from Holder's inequality 

and the continuity of the trace operator from iJ5 +e (f2) to L2(T). 

The exact solution for each example can be expressed by an infinite series involv- 
ing the Mittag-Leffler function E a: i(z). To accurately evaluate the Mittag-Lcfilcr 
functions, we employ the algorithm developed in |13j . To discretize the problem, we 
divide the unit interval (0, 1) into N — 2 k equally spaced subintervals, with a mesh 
size h — 1/N so that [0, l] 2 is divided into N 2 small squares. We get a symmetric 
mesh for the domain [0, l] 2 by connecting the diagonal of each small square. All 
the meshes we have used are symmetric and therefore both semidiscrete Galerkin 
FEM and lumped mass FEM have the same theoretical accuracy. Unless otherwise 
specified, we have used the lumped mass method. 

To compute a reference (replacement of the exact) solution we have used two 
different numerical techniques on very fine meshes. The first is based on the exact 
representation of the semidiscrete lumped mass solution Uh by 

N-l 

U h (t)= ]r £a,l(-A£ >m nK< m )< m , 
n,m=l 

where m (x, y) = 2 sin(n7ra;) sin(m7ry), n, m = 1, . . . , N — 1, with x, y being grid 
points, are the discrete eigenfunctions and 

h 4 / 2 nnh 2 Tn'^h\ 

X n, m = I Sin ~Y~ + sm ~Y~ ) 

are the corresponding eigenvalues. 

The second numerical technique is based on fully discrete scheme, i.e., discretiz- 
ing the time interval [0, T] into t n = nr, n = 0, 1, . . . , with r being the time step 
size, and then approximating the fractional derivative d^u(x, t n ) by finite difference 
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m 

t-n—j ) t n —j — i J 



(5.1) d?u{x,t n )^ T{2 a) J2 b r 

where the weights bj = (j + 1) 1_Q — j 1 ^", j = 0, 1, . . . ,n — 1. This fully discrete 
solution is denoted by Uh- Throughout, we have set r = 10~ 6 so that the error 
incurred by temporal discretization is negligible (see Table [6] for an illustration) . 

We measure the accuracy of the approximation Uh(t) by the normalized error 
\\u(t) — Uh(t)\\/\\v\\ and ||V(u(f) — Ufc(t))||/||u||- The normalization enables us to 
observe the behavior of the error with respect to time in case of nonsmooth initial 
data. 

Smooth initial data: example (a). In Table [4] we show the numerical results for 
t = 0.1 and a — 0.1, 0.5, 0.9. Here ratio refers to the ratio between the errors as 
the mesh size h is halved. In Figure [IJ we plot the results from Table [4] in a log-log 
scale. The slopes of the error curves are 2 and 1, respectively, for L^- and 7J 1 -norm 
of the error. This confirms the theoretical result from [B]. 



Table 4. Numerical results for smooth initial data, example (a), 
t = 0.1. 



a 


h 


1/8 


1/16 


1/32 


1/64 


1/128 


ratio 


rate 


0.1 


Z/2-norm 


9.25e-4 


2.44e-4 


6.25e-5 


1.56e-5 


3.85e-6 


« 4.01 


0(h 2 ) 


i/ i -norm 


3.27c-2 


1.66c-2 


8.40e-3 


4.21c-3 


2.11e-3 


ps 1.99 


0(h) 


0.5 


i/2-norm 


1.45e-3 


3.84e-4 


9.78e-5 


2.41e-5 


5.93e-6 


PS 4.02 


om 


if x -norm 


5.17e-2 


2.64e-2 


1.33e-2 


6.67e-3 


3.33e-3 


ps 1.99 


0(h) 


0.9 


L 2 " n orm 


1.88e-3 


4.53e-4 


1.13e-4 


2.82e-5 


7.06e-6 


ps 4.00 


0(fc a ) 




6.79e-2 


3.43e-2 


1.73e-2 


8.63e-3 


4.31e-3 


ps 2.00 


0(h) 




Figure 1. Error plots for smooth initial data, Example (a): a — 
0.1,0.5,0.9 at t = 0.1. 



Intermediate smooth data: example (b). In this example the initial data v(x) is 
in Hq(£1) and the numerical results are shown in Table [5] The slopes of the error 
curves in a log-log scale are 2 and 1 respectively for L 2 - and iJ 1 -norm of the errors, 
which agrees well with the theory for the intermediate case [BJ. 



GALERKIN FEM FOR FRACTIONAL PDE'S WITH NON-SMOOTH DATA 

Table 5. Intermediate case (b) with a = 0.5 at t = 0.1. 



li 



h 


1/8 


1/16 


1/32 


1/64 


1/128 


ratio 


rate 


L2-error 


3.04c-3 


8.20e-4 


2.12e-4 


5.35e-5 


1.32e-5 


w 3.97 


0(h z ) 


i? 1 -error 


5.91e-2 


3.09e-2 


1.56e-2 


7.88e-3 


3.93e-3 


w 1.98 


0(h) 



Nonsmooth initial data: example (c). First in Ta ble [6| we compare fully discrete so- 
lution Uh via the finite difference approximation (5.1 ) with the semidiscrete lumped 
mass solution Uh via eigenexpansion to study the error incurred by time discretiza- 
tion. We observe that for each fixed spatial mesh size h, the difference between 
Uh, the lumped mass FEM solution, and Uh decreases with the decrease of r. In 
particular, for time step r = 10 -6 the error incurred by the time discretization is 
negligible, so the fully discrete solutions Uh could well be used as reference solutions. 
In Table [7] and Figure [2] we present the numerical results for problem (c). These nu- 



Table 6. The difference Uh — Uh, nonsmooth initial data, example 
(c): a = 0.5, t = 0.1 



Time step 


h 


1/8 


1/16 


1/32 


1/64 


1/128 


T = 10- 2 


Z,2-norm 


2.03e-3 


2.016-3 


2.00e-3 


2.00e-3 


2.00e-3 


i? i -norm 


9.45e-3 


9.17e-3 


9.10c-3 


9.08e-3 


9.07e-3 


T = 10~ 4 


L2-norm 


1.81e-5 


1.79e-5 


1.79e-5 


1.79e-5 


1.796-5 


i? 1 -norm 


8.47e-5 


8.22e-5 


8.15e-5 


8.13e-5 


8.13e-5 


T = 10" 8 


Z/2-norm 


1.80e-7 


1.786-7 


1.78&-7 


1.78&-7 


1.78e-7 


i? 1 -norm 


8.42e-7 


8.17e-7 


8.10e-7 


8.10e-7 


8.09e-7 



merical results fully confirm the theoretically predicted rates for nonsmooth initial 
data. 











jay /, p 








2 ^ 




/ / / 
/ / / , 
$ a sf / 
/ / / / 
/ / # / i 

V / 

/ / / 1 — 1 — 

* <2> 




10- 2 




1 y 







10 



10" 



10 



-*-H',t=0.001 



- & -H',t=0.01 



"H '.1=0.1 



Figure 2 . Error plots for lumped FEM for nonsmooth initial data, 
Example (c): a — 0.5. 



Very weak data: example (d). The empirical convergence rate for the weak data Sr 
agrees well with the theoretically predicted convergence rate in Theorem 3.1 which 
gives a ratio of 2.82 and 1.41, respectively, for the L2- and i? 1 -norm of the error; 
see Table [9] Interestingly, for the standard Galerkin scheme, the L 2 -norm of the 
error exhibits super-convergence; see Table [8j 
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Table 7. Error for the lumped FEM for nonsmooth initial data, 
example (c): a = 0.5 



Time 


h 


1/8 


1/16 


1/32 


1/64 


1/128 


ratio 


rate 


t = 0.001 


Z/2-norm 


1.55e-2 


3.99e-3 


1.00e-3 


2.52e-4 


6.26e-5 


ps 4.01 




i/ 1 -norm 


6.05e-l 


3.05e-l 


1.48e-l 


7.29e-2 


3.61e-2 


ps 2.00 


0(A) 


t = 0.01 


L 2 -norm 


8.27e-3 


2.10e-3 


5.28e-4 


1.32e-4 


3.29c-5 


ps 4.01 




H i -norm 


3.32e-l 


1.61e-l 


7.90e-2 


3.90e-2 


1.93e-2 


ps 2.02 


o(M 


t = 0.1 


L2-norm 


2.12e-3 


5.36e-4 


1.34e-4 


3.36e-5 


8.43e-6 


ps 3.99 




H 1 -norm 


8.23c-2 


4.01e-2 


1.96e-2 


9.72c-3 


4.84e-3 


ps 2.01 


0(/i) 



Table 8. Error for standard FEM: initial data Dirac 8- function, 
a = 0.5 



Time 


h 


1/8 


1/16 


1/32 


1/64 


1/128 


ratio 


rate 


t = 0.001 


L2-norm 


5.37c-2 


1.56e-2 


4.40e-3 


1.23c-3 


3.41e-4 


« 3.57 


0(0**) 


H -norm 


2.68e0 


1.76e0 


1.20e0 


8.21e-l 


5.68e-l 


ps 1.45 


0(0) 


t = 0.01 


L2-norm 


2.26c-2 


6.20e-3 


1.67e-3 


4.46e-4 


1.19e-4 


ps 3.74 


0(h LM ) 


H -norm 


9.36e-l 


5.90e-l 


3.92e-l 


2.65e-l 


1.84e-l 


ps 1.46 


0(0) 


t = 0.1 


L2-norm 


8.33c-3 


2.23e-3 


5.90e-3 


1.55c-3 


4.10e-4 


ps 3.77 


O(0 An ) 


H -norm 


3.08e-l 


1.91e-l 


1.26e-l 


8.44c-2 


5.83c-2 


ps 1.46 


0(0) 



Table 9. Error for lumped mass FEM: initial data Dirac 8- 
function, a = 0.5 



Time 


h 


1/8 


1/16 


1/32 


1/64 


1/128 


ratio 


rate 


t = 0.001 


L2-norm 


1.98e-l 


7.95e-2 


3.00e-2 


1.09e-2 


3.95e-3 


ps 2.75 


0(0) 


H -norm 


5.56e0 


4.06e0 


2.83e0 


2.02e0 


1.41e0 


ps 1.42 


0(0) 


t = 0.01 


L2-norm 


6.61c-2 


2.56e-2 


9.51e-3 


3.47e-3 


1.25c-3 


ps 2.78 


0(0) 


if 1 -norm 


1.84e0 


1.30e0 


9.10e-l 


6.40e-l 


4.47e-l 


ps 1.42 


0(0) 


t = 0.1 


Z/2-norm 


2.15e-2 


8.13e-3 


3.01e-3 


1.09e-3 


3.95e-4 


ps 2.75 


0(0) 


H -norm 


5.87e-l 


4.14e-l 


2.88e-l 


2.03c-l 


1.41e-l 


« 1.43 


0(0) 













" t f f 






t=0.001 




; 1 


Lb y 






i i i 

4 a <t 






-»-H 1 ,t=0.001 












9 9 9/ 


2 




-e-L t=0.01 


10 2 










/ / / / 
d <f 




- B "H 1 ,t=0.01 
-o-L 2 ,1=a.1 


10" 4 


10 3 


10 2 


10" 1 


10° 10 1 




-»-H 1 ,t=0.1 



error 



Figure 3. Error plots for Example (d): initial data Dirac 8- 
function, a = 0.5. 
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